//==============================================================================
// LSLAligner.cc
// 
// Linear Space Local Alignmnet. 
// 
//==============================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8
//

#include "LSLAligner.hh"

//long limitingCaseCount = 0;


//------------------------------------------------------
// * GetAlignment
//------------------------------------------------------
// Recovers the best locally aligned subsequence of two 
// sequences.  
//
void
LSLAligner::GetAlignment(Sequence &Source,Sequence &Query)
{
  mCols = Source.size()+1;
  mRows = Query.size()+1;

  mS = Source;
  mQ = Query;

  // Allocate two columns for score table and two columns 
  // for intersect forwarding table ("c" table). 
  mScoreTable.Allocate(2,mRows); 
  mCTable.Allocate(2,mCols);

  // Make sure it is initialized to something sane
  for(int r = 0;r < 2;r++){
	 for(int c = 0; c < mCols;c++){
		mScoreTable[r][c] = 0;
	 }
  }

  for(int r = 0;r < 2;r++){
	 for(int c = 0; c < mCols;c++){
		mCTable[r][c] = 0;
	 }
  }

  mMaxScore = -99999;
  mMaxRow = 0;
  mMaxCol = 0;

  //  cout<<"GetAlignment(00RC) "<<0<<"\t"<<0<<"\t"<<Source.size()-1<<"\t"<<Query.size()-1<<endl;
  GetAlignment(1,1,Query.size(),Source.size());
}

//------------------------------------------------------
// * GetAlignment
//------------------------------------------------------
// Divide and conquer approach to alignment.  Finds a 
// single point on the optimum alignment (u,v) and
// recursively splits the alignment problem around that
// point.  When the sub-problem reaches the scale of aligning
// a single pair of letters (or letters and gaps), the 
// alignment is inferred directly from the shape of the 
// sub-problem (2x1, 1x2, or 2x2, for gap-letter, letter-gap, 
// and letter-letter respectively).  The full alignment is 
// simply the concatenation of the alignments from the left
// and right half of each sub-problem.  Here, I construct 
// the alignment by right-end concatenation on the mAlignS, and
// mAlignQ vectors, which works out because recursion is depth first 
// on the left half then right.  
//
void 
LSLAligner::GetAlignment(int r1,int c1,int r2,int c2)
{

  cout<<"GetAlignment(r1c1r2c2)\t"<<r1<<"\t"<<c1<<"\t"<<r2<<"\t"<<c2;

  int w = c2-c1+1;
  int h = r2-r1+1;
  int u,v;
  char ans;

  cout<<"\tw="<<w<<"\th="<<h<<endl;

  // Sanity check
  if ((w<0) || (h<0)){
	 cout<<"Invalid range.\n";
	 exit(1);
  }

  // 1x1, 2x1, 1x2 are all limiting cases.  2x2 is not. 
  //if (((w <= 2) && (h <= 2)) && !(w==2 && h ==2)){
  if ((w < 2) || (h < 2)){
	 HandleLimitingCase(r1,c1,r2,c2);
  }else{
	 // The midpoint row.  
	 u = ((r2-r1)/2)+r1; 
	 // Find the corresponding column where the best path 
	 // intersects this row. 	
	 cout<<"Find Intercept (u,r1,c1,r2,c2): "<<u<<"\t"<<r1<<"\t"<<c1<<"\t"<<r2<<"\t"<<c2<<endl;
	 v = FindInterceptWithBestPath(u,r1,c1,r2,c2);

	 cout<<"LEFT: ";
	 cout<<"u="<<u<<"\tv="<<v<<"\t";
	 GetAlignment(r1,c1,u,v);

	 cout<<"RIGHT: ";
	 cout<<"u="<<u<<"\tv="<<v<<"\t";
	 GetAlignment(u+1,v+1,r2,c2);
  }
}

//------------------------------------------------------
// * HandleLimitingCase
//------------------------------------------------------
// Using push_back only works if recursion solves left 
// hand side before right hand side.
//
void
LSLAligner::HandleLimitingCase(int r1,int c1,int r2,int c2)
{

  int w = c2-c1+1;
  int h = r2-r1+1;

  cout<<"Handle limit "<<r1<<"\t"<<c1<<"\t"<<r2<<"\t"<<c2<<" w="<<w<<" h="<<h<<endl;

  if ((w == 1) && (h == 1)){
	 mAlignS.push_back(mS[c2-1]);
	 mAlignQ.push_back(mQ[r2-1]);

	 cout<<mAlignS<<endl;
	 cout<<mAlignQ<<endl;
	 return;
  }

  if ((w==1) && (h > 1)){
	 for(int i = h-1;i >= 0;i--){
		mAlignS.push_back('-');
		mAlignQ.push_back(mQ[r2-1-i]);
	 }
	 cout<<mAlignS<<endl;
	 cout<<mAlignQ<<endl;
	 return;
  }


  if ((w > 1) && (h == 1)){
	 for(int i = w-1;i>= 0;i--){
		mAlignS.push_back(mS[c2-1-i]);
		mAlignQ.push_back('-');
	 }
	 cout<<mAlignS<<endl;
	 cout<<mAlignQ<<endl;
	 return;
  }


}

/*
Old handle limit:

 if (((w==1) && (h == 1)) ||
      ((w==0) && (h == 0))){
	 //cout<<"Case 1"<<endl;
	 mAlignS.push_back(mS[c2-1]);
	 mAlignQ.push_back(mQ[r2-1]);
  }else{
	 if ((w== 0) && (h == 1)){
		//cout<<"Case 3"<<endl;
		mAlignS.push_back('-');
		mAlignQ.push_back(mQ[r2-1]);
	 }else{
		if ((w==1) && (h == 0)){
		  //cout<<"Case 2"<<endl;
		  mAlignS.push_back(mS[c2-1]);
		  mAlignQ.push_back('-');
		}else{
		  cout<<"Invalid limiting case detected."<<endl;
		  cout<<"w="<<w<<" h="<<h<<"\n";
		  exit(1);
		}
	 }
  }


 */


//------------------------------------------------------
// * FindInterceptWithBestPath
//------------------------------------------------------
// Finds the intersection of the optimum path with column
// u.  This uses the "c" table method described in Durbin, Eddy, 
// Krough, and Mitchison, "Biological Sequence Analysis".
//
int
LSLAligner::FindInterceptWithBestPath(int u,int r1,int c1,int r2,int c2)
{
  
  int diagScore,leftScore,topScore,temp;
  int rows = r2-r1+1;
  int cols = c2-c1+1;
  int current = 1,previous = 0;   // score table columns
  int i,j;
  
  // Each time we do this we start as a virgin. 
  mMaxRow = 0;
  mMaxCol = 0;
  mMaxScore = -999;
  
  // Adjust u to local sub-table coordinates 
  u = u - r1+1;
  
  // Fill in initial column with gap scores and the C-table with 
  // row indices as initial lables. 
  for(int c = 0;c<=cols;c++){
	 mScoreTable[previous][c] = 0;
	 mScoreTable[current][c] = 0;
  }

  for(int c = 0;c<=cols;c++){
	 mCTable[previous][c] = 0;
	 mCTable[current][c] = 0;
  }  

  int maxRowIdx = 0,maxColIdx = 0;

  // Compute optimum alignment score one row at a time for 
  // the given range.  
  int r,c,rIdx=0,cIdx = 1;
  do{
	 rIdx++;
	 
	 // Fill in the beginning of each row with gap scores.
	 mScoreTable[current][0] = 0;
	 
	 for(cIdx = 1;cIdx <= cols;cIdx++){
		
		//		cout<<"Score at (r,c): "<<rIdx<<","<<cIdx<<"\tSource["<<c1+cIdx-2<<"]=";
		//cout<<mS[c1+cIdx-2]<<"\t";
		//cout<<" Query["<<r1+rIdx-2<<"]=";
		//cout<<mQ[r1+rIdx-2]<<"\t ==>";
		int  matchScore = Score((char)mS[c1+cIdx-2],mQ[r1+rIdx-2]);
		//cout<<matchScore<<endl;
		
		diagScore = mScoreTable[previous][cIdx-1] + matchScore;
		topScore = mScoreTable[previous][cIdx] + mGapScore;
		leftScore =  mScoreTable[current][cIdx-1] + mGapScore;	
		temp=max(diagScore,leftScore);
		temp=max(temp,topScore);		
				
		//		mScoreTable[rIdx][current] = temp;
		
		// With local alignment, we can't assume that the best overall score
		// will end up in the corner, it might be anywhere. 
		if (temp > mMaxScore){
		  mMaxScore = temp;
		  maxRowIdx = rIdx;
		  maxColIdx = cIdx;
		}
		
		//		cout<<"left="<<leftScore<<"\ttop="<<topScore<<"\tdiag="<<diagScore<<endl;
		if (temp >=0){
		  if (rIdx == u){
			 if (temp == diagScore) c =cIdx-1; 			
			 if (temp == topScore) c = cIdx;
			 if (temp == leftScore) c = cIdx -1;
		  }else{
			 if (temp == diagScore) c = mCTable[previous][cIdx-1];
			 if (temp == topScore) c = mCTable[previous][cIdx];
			 if (temp == leftScore) c = mCTable[current][cIdx-1];
		  }
		  mCTable[current][cIdx]=c;
		  mScoreTable[current][cIdx]= temp;
		}else{
		  mCTable[current][cIdx]=0;
		  mScoreTable[current][cIdx]= 0;
		}
	 }		  

	 //	 cout<<endl;
	 //PrintCTable(cols,previous,current);
	 //PrintScoreTable(cols,previous,current);

	 Swap(previous,current);
	 
	 if (!(rIdx%1000)) cout<<rIdx<<endl;
	 
  }while(rIdx < rows);
  //  while((rIdx + r1) <=r2);  	


  // Undo last,unused swap. 
  Swap(previous,current);  


  cout<<endl;
  PrintCTable(cols,previous,current);
  PrintScoreTable(cols,previous,current);

  int rValue = c1+mCTable[current][maxColIdx];
  cout<<"c1="<<c1<<" maxColIdx="<<maxColIdx<<"\trValue="<<rValue<<endl;
  cout<<"========================================================================"<<endl;
  cout<<"========================================================================"<<endl;

  return(rValue);
}

void 
LSLAligner::PrintCTable(int cols,int previous,int current){
  cout<<"CTable:"<<endl;
  for(int c = 1; c<= cols;c++){
	 cout<<mCTable[previous][c]<<" ";
  }
  cout<<endl;
  for(int c = 1; c<= cols;c++){
	 cout<<mCTable[current][c]<<" ";
  }
}

void 
LSLAligner::PrintScoreTable(int cols,int previous,int current){
  cout<<endl;
  cout<<"ScoreTable:"<<endl;
  for(int c = 1; c<=cols;c++){
	 cout<<mScoreTable[previous][c]<<" ";
  }
  cout<<endl;
  for(int c = 1; c<= cols;c++){
	 cout<<mScoreTable[current][c]<<" ";
  }
  cout<<endl;
}



